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Abstract 

The muon intensity attenuation method to detect heterogeneities in large matter 
volumes is analyzed. Approximate analytical expressions to estimate the collection 
time and the signal to noise ratio, are proposed and validated by Monte Carlo 
simulations. Important parameters, including point spread function and coordinate 
reconstruction uncertainty are also estimated using Monte Carlo simulations. 
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1 Introduction 



The first practical application of cosmic muon intensity attenuation measure- 
ments dates back 60 years, when attempts were made to evaluate snow layers 
in Australian mountains [1]. Ten years later Alvarez et al. [2] used this princi- 
ple to search for hidden chambers in the great Egyptian pyramid of Cheops, in 
Giza. The use of near-horizontal muon attenuation to probe active volcanoes 
[3-5] or to detect and locate heavy metals in border-crossing vehicles using 
scattering muon tomography [6, 7], which can also be combined with muon di- 
rectional intensity attenuation, electro-magnetic shower induction, and muonic 
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X-ray emission as additional means for the non-destructive assay of small and 
medium size objects [8], or the use of muon-induced neutron emission to gauge 
carbon liner thickness in steel industry blast furnaces [9], are good examples 
of the continued interest in applying cosmic muon detection to a wide variety 
of practical problems. 

In our previous works [10,11], we explored the viability of carrying out an 
experiment similar to that of Alvarez et al.[2] in the Pyramid of the Sun at 
Teotihuacan, in Central Mexico. Here we evaluate important experimental and 
simulation planning parameters, such as detector size, coordinate resolution 
and point spread function (PSF), the signal-to-noise ratio and observation time 
to detect a structure of a given size. Some of those parameters are estimated 
analytically using an approximate description of muon directional intensity 
(DI). Here we shall concentrate on an energy range relevant to explore matter 
depths in the range 5 — 50 kg /cm 2 . Monte Carlo simulations have also been 
performed to estimate the magnitude of PSF and the coordinate resolution as 
functions of the matter thickness along the muon track, as well as to test the 
approximate expression proposed below for the signal to noise ratio. 



2 Muon energy spectrum and signal to noise ratio 

The muon DI attenuation method [1] requires a knowledge of the initial cosmic 
muon differential energy spectrum at the site of interest F(E, 6), where E is 
the muon energy and G the zenith angle. For the vertical direction (0 = 0) 
at sea level F can be represented by [12]: 

F(E,Q = 0) = GE- n ^[cm- 2 s- 1 sr~ 1 GeV- 1 ], (1) 

where G, is a normalization constant and the spectral index n(E) is a slowly 
varying energy- dependent function. A well established[13,14] reference value 
for the spectral index is n[E = 100 GeV) = 2.7. This fact can be used to 
extract the value of G at that point. The energy dependence of n(E) can then 
be determined using the experimental data. 

Muon differential energy spectrum measurements exist with sufficient accuracy 
at sea level and within the energy range E < 100 GeV [15]. These data have 
been compiled and analyzed by various authors[13-16], resulting in a variety 
of parameterizations. Among them the one by Hebbeker and Timerman [15] 
describes the data in the energy region of 10 — 100 GeV with an accuracy 
better than 2%, while for E > 100 GeV the experimental data has larger 
uncertainties, of the order of 10% — 15%. The results of n(E) obtained using 
the parameterizations [15] and [13,14] is represented in Fig. 1. 
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The DI attenuation method [1] is based on using the muon integral energy 
spectrum <&(E > El, 0) (where El is the muon mean threshold energy to 
cross the volume of interest expressed in GeV, for a matter thickness L along 
the muon track, given in density length units [3]), defined as [12]: 



$>(E L , = 0) = J F(E, = 0)dE. (2) 

E L 



As shown in Fig 1 (solid line), n(E) varies less than ±1.5% in the energy 
region E = 10 — 100 GeV, so that an n(E) = n c = constant approximation 
can be used to obtain a simplified analytic expressions for experiment planning 
purposes. Using n c = 2.687 provides (see Fig 2) a ~ 10% accuracy for both 
differential and integral vertical energy spectra. 

Locating an heterogeneity in the matter volume using muon intensity atten- 
uation is based on calculating the difference AiV between the experimental 
number of muons iV (after crossing the volume of interest) and the simulated 
number of muons in an equivalent volume in which no heterogeneity is in- 
cluded. Hereafter we shall refer to this as the " simulated background" . If o"av 
represents the uncertainty of AiV, the signal to noise ratio, (£) is defined as: 

( ^ ^ (3) 

Note that the uncertainty associated with the simulated background can all- 
ways be reduced below the experimental one, thus a^N ~ VN- 

Let us address the question of what is best to measure, F or Within n(E) = 
n c approximation, one can calculate the time-integrated vertical muon count 
iV and the difference AiV for an energy E L , and the energy difference AE 
associated to the heterogeneity, registered in a detector having surface S, 
solid angle acceptance Q and during collection time T using iV = GQ(El, = 
0)ttST and AiV ps GF(E l , = 0)AEilST (assuming that AE < E L ). Thus, 
Eq.(3) can be used to estimate the signal to noise ratio associated with 
the integral intensity $. A similar procedure can be used to calculate the signal 
to noise ratio (£p) associated with a finite energy interval e <C E L using F. 
Then the ratio between and £p can be written approximately as: 



\ 



E L (n c -l) ^ 



en\ 



From this expression it is clear that it becomes preferable to measure the 
integral intensity (rather than the differential one) when E L > en 2 c / (n c — 1) w 
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4.28e. Since the n(E) = n c approximation is valid for E > lOGeV, considering 
£ = l GeV (to fulfill e < E L ) implies E L > 4.28 GeV for the > £ F 
condition to be fulfilled. Still, the use of differential energy spectrum gives a 
better ratio for AN/N [2]. This consequence of the power-law shape of the 
muon energy spectra justifies the important experimental simplification of 
measuring $ rather than F. Thus, from now on we consider only the signal 
to noise ratio for the integral energy intensity. 

Taking into account that there is a unique relationship between a moun energy 
and the mean range, £ can also be expressed as a function of L, which is 
more convenient in some practical applications. For standard soil matter, with 
average density 2 x 10 -3 kg /cm 3 , and within the energy interval 5 — 100 GeV 
the calculated data for the range [17] can be fitted using E L = cL k 1 , where L 
is expressed in kg /cm 2 , c = 1.84 GeV / (kg / cm 2 ) k and k = 1.074. 

When the localized density heterogeneity has a size AL <C L, and still within 
the n(E) = n c approximation £ can be expressed as a function of AL and L 
as follows: 

^^ r "(^- 1 ) c( n e -l)/2^ e -l)/2 + l . (5) 



This implies that, for a fixed heterogeneity size AL, £ decreases as L~ im& . To 
illustrate the impact of our approximation on £, the ratio £, apP /£,ex between £ 
corresponding to the approximate (£ apP ) and the exact (£ ex ) differential muon 
intensities, for different AL/ L values, as a function of L is presented in Fig 3. 
There is an agreement within 10% for AL/ L < 10% in the range L = 8 — 50 
kg /cm 2 . 

Concerning the experimental statistics necessary to achieve a given £, still 
within the constant spectral index n(E) approximation, we find: 



Thus, the required collection time for a given £ can be estimated using: 



t1 n n c -\ Tk(n c -l)+2 

2 (7) 

~ GnS(n c -l)k 2 AL 2 ' 1 ; 

1 This is an overestimation of El, because the steepness of F oc E~ 2 ' 7 leads to 
an asymetric distortion, an enhancement of the low energy contribution. Hear we 
neglect this effect due to the narow energy distribution of muons having a given 
range value L 
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Since the collection time increases as L 3,812 , to minimize T it is necessary to 
choose a detector location where L gets its minimum possible value for the 
region of interest. 

Using the value of AL, on can estimate the geometrical length AX of the 
heterogeneity, measured along the muon direction as: 

AX = AL/p, (8) 



where p is the average matter density, expressed in kg /cm 3 . Assuming that 
the density heterogeneity length AX <C X, where X is the total length of the 
muon trajectory inside the volume of interest, the uncertainty of AX, cr^x 
can be estimated from: 

« + Cf) 2 + ((kn c y + (k- im^y + (i + (kn c y)(^y,(9) 



where ax and a p are the uncertainties corresponding to the total length X 
and mean density p, respectively. Note that the first term in Eq.(9) depends 
only on the statistics. In the next section we estimate the other terms in this 
equation. 

To apply these equations for 0^0, the differential energy intensity for a 
given angle should also be parameterized in the form given by Eq. (1). Here, to 
describe the angular variation we used cos(Q) m ( E \ so the differential intensity 
for a given becomes: 

F(E, 0) = GE~ n W cos(0) m(£) , (10) 



Because of the poor quality of the experimental data at large angles, m(E) 
was obtained by fitting CORSIKA[18] simulated data. The procedure was 
to simulate the data over the entire angular range and normalize the results 
to the data in the small angular region < 5°. This provided an educated 
guess for the polynomial fit to the large angular part. The obtained results of 
m(E) at sea level were fitted by a polynomial function, resulting in m(E) = 
1.53 - 0.0484E + 0.000545E 2 - 2.79 x 10~ 6 E 3 + 5.11 x 1(T 9 £ 4 , where E is 
expressed in GeV. This is represented in Fig 4 (solid line) together with the 
experimental data. As can be seen, the obtained results agree with the data 
within the experimental errors. 

To present Eq (10) in the form of Eq (1) we introduced a function g(Q) = 
Gcos(0) m ^ =loGey ) to describe the angular dependence of G, and rj(E, 0) for 
the spectral index, which becomes n(E) at = 0. 

F(E,Q) = g(Q)E-^ E ' e \ (11) 
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For a given value of 0, approximate analytic calculations can also be carried 
out using a constant spectral index 1] C (Q), which likewise becomes n c at = 0. 
The value of rj c (Q) would depend on the zenith angle, because of the additional 
cos(6) m( ^ dependence. 

As shown in Fig 2, within the angular range = 0° — 60°, i] c (Q) varies between 
2.687 and 2.53, respectively. Also, within a 10% accuracy, at O = 60° the 
approximation is appropriate for energies larger than 25 GeV. The decrease 
in r] c (Q) relative to the vertical direction can reach 5 — 6% and almost twice 
in g(Q) parameter for = 60°. For ?? c (0) a simple parameterization rj c (Q) = 
n c cos(#) ' 0907 can be used. Another problem is that most of the available data 
compilations correspond to sea level measurements, while for other altitudes 
the experimental data is very scarce. For small altitude changes (h < 1000 
m, where h in m measured from the sea level), the variation in the shape of 
vertical differential energy spectrum for > 10 GeV can be neglected. This 
is because the contribution of muon decay for those energies is negligible [3]. 
Thus the altitude correction to the differential intensity can be applied for the 
normalization coefficient in Eq (1), using the exponential expression proposed 
in Ref. [15,19] as: 

G h = Gexp(h/h (E L = IQGeV)), (12) 

Where Gh is the normalization coefficient for the altitude h and ho is a function 
of El defined in [15]. For large altitudes (for example h = 2300 m), the m(E) 
obtained by fitting CORSIKA simulation data is shown in Fig 4 (dashed line). 
From the figure one can observe a difference between the m(E) values obtained 
for the sea and Mexico City levels. So, for large altitudes, the value of rj c (Q) 
and g(Q) should be corrected. 



3 Uncertainties 

The second term in Eq.(9) is the uncertainty of the muon differential energy 
spectrum F. The uncertainty associated with the vertical differential energy 
spectrum has already been discussed in Section 2. Concerning cos(0) m ^, 
Eq.(10) provides good results at sea level, better than 5%, within 0° < < 
40°. As illustrated in Fig. 5, beyond that angular region there is less data, and 
the corresponding error bars are bigger. For larger angles, the agreement with 
the experiment is ~ 10% — 20% for the energy interval E < 100 GeV (Fig. 6). 

Thus, for non-sea-level applications we also relied on CORSIKA simulations 
to estimate the altitude behavior of the muon intensity. For altitudes h < 2300 
m the accuracy of Eq. (12) has been checked using simulation data resulting 
in values better than < 5 %. 
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The third term in Eq.(9) represents the geometrical uncertainties, associated 
with the description of the external shape of the investigated volume, as well 
as those associated with the detector location. They define the muon path 
length X, which is known with some uncertainty ax- The AL determination 
is affected by these geometrical uncertainties. However, the location of the 
detector only contributes to the systematic errors, introducing a fixed dis- 
placement but with no fluctuating character. 

The last term in Eq.(9) represents the uncertainty in the internal density dis- 
tribution. The simulation requires the best possible knowledge of the density 
and chemical composition of the materials contained in the volume of interest. 
This can be achieved by sampling the volume along some significant direction. 
When this procedure is carried on systematically, one can define a sampling 
length X samp , and estimate the density uncertainty a s p amp . Yet, we don't ex- 
pect randomly distributed density variations to be important, because over a 
long X the uncertainty in the mean value of the density depends on X, i.e., 



Finally, there are a number of smaller factors known to affect the cosmic muon 
intensity (E > 10 GeV), such as the Solar modulation, the Earth position 
and even the local atmospheric conditions, temperature and pressure. The 
importance of these effects has been reported by some authors [15,12,20] to 
be of the order 1%. 



4 Limitations 

With the uncertainties described in the previous section, now we can estimate 
the minimum heterogeneity size as the uncertainty of X (assuming negligible 
statistical error), using: 



One factor not yet included in these estimations is the possibility that the 
heterogeneity of interest may not be a simple localized density excess, or defect, 
relative to the mean density calculated along a given muon trajectory, but a 
combination of them, reducing the sensitivity of the method. To quantify this, 
let us define the AX X and AX 2 as the total lengths corresponding to the 
densities p\ and p 2 , respectively, where p\ < p < p 2 . Taking into account that 
£ oc AL = p 1 AX 1 + p 2 AX 2 — p(AX 1 + AX 2 ) Eq. (5), one can conclude that 
£ is reduced when p approaches (p\AXi + p 2 AX 2 ) j '(AX 1 + AX 2 ) 




AX,, 



mm 



(X/k) y/(kc x /X)* + (* P (X)/P) 2 + (°f/F)\ 



(13) 
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5 PSF and coordinate reconstruction uncertainty 

The effects of muon multiple scattering are best studied using Monte Carlo 
simulations. For this purpose we use the simulation package GEANT4 [24] 
which allows us to closely reproduce the relevant physical processes. An im- 
portant parameter for heterogeneity localization is the coordinate reconstruc- 
tion uncertainty (ay), associated with the reconstructed muons location at a 
given density length distance from the detector (L d ). The knowledge of the 
coordinate transform (displacement) uncertainty (a p ) is important for the op- 
timization of experimental equipment. For example, if the detector size is less 
than 4o" p , then important information for image reconstruction will be lost. 
The relation between the object lateral size and a p is also important. If the 
object's lateral size is comparable with a p , then one should expect a decrease 
in f. 

The above-mentioned parameters depend not only on Ld but also on the total 
thickness of the absorber for the given direction, because of the shift in the 
corresponding energy region, resulting a in a reduction of the fraction of low 
energy muons having large multiple scattering. 

The Ld dependence of the uncertainties oy and a p will be estimated within 
the range 2 — lQkg/cm 2 and for a total density length of 16kg/cm 2 , which is 
expected to be the maximum thicknesses of our trial experiment [10] (see next 
section). To estimate o r and a p , a vertical beam of muons with the energy 
distribution given by Eq (1) was simulated, passing through different matter 
thicknesses, with the mean density p ~ 2 x 10~ 3 kg/cm 3 [10] and the standard 
soil material [23]. The penetrated muons were detected by the detector located 
below the matter volume having an ideal resolution. 

The value of a p is estimated as the standard deviation of the distribution of 
coordinate displacements on the detector surface. The back-projected coordi- 
nate distribution on the volume surface is obtained using the detected track 
information and a r is estimated as the standard deviation of the corresponding 
distribution. 

The results of a r and a p are presented in Fig. 7. Hence mentioned by Alvarez et 
al. [2], the use of low energy cuts on the detected muons does not significantly 
improve the resolution parameters. Both parameters have approximately linear 
dependence on the density length. In the figure the parameters for the case 
when L d = L are also presented, which clearly demonstrate the effect of total 
thickness. 

In Ref. [25] we showed that these simulations for multiple scattering angle 
with large energy cuts, and for different density lengths, are in agreement 
with the di-muon angular distribution in underground measurements by the 
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L3C Collaboration [26]. 



6 Heterogeneity simulation 

The purpose now is to validate the obtained approximate relations using re- 
alistic muon energy distributions based on Eq (10), and including effects of 
multiple scattering. To illustrate this we use GEANT4, taking the Pyramid 
of the Sun in Teotihuacan, Mexico, as a model of the matter volume. The 
complex geometrical shape of this monument is shown in Fig 8. 

Two simulations using Eq. (10) have been performed with, and without, in- 
ternal structures in the volume of the pyramid. For each simulation 5 x 10 7 
events for zenith angle — 50 were generated. As model heterogeneities we 
considered two empty volumes (see Fig 8). A chamber having xyz dimensions 
(16 x 5 x 4m 3 ), which are much larger than a p , so that multiple scattering 
will not strongly affect its detection. The other heterogeneity considered is a 
60 x 2 x lm 3 tunnel, which is twice longer than an existing modern structure 
known as the Smith Tunnel. In this case the width of the tunnel is approxi- 
mately the triple of a p and one should expect a lower value for £. The detector 
sensitive area of 1 x lm 2 is also slightly smaller than 3a p , and will introduce 
some decrease of £. 

The results of the Monte Carlo simulations for the detected muons, with hypo- 
thetical empty volumes are presented in Fig. 9 showing a plot of the projection 
angles Q x and O z where x and z are the coordinates in the horizontal plane. 
The chamber is clearly seen, while a hint of the tunnel can also be observed. 
To make the latter more visible, the £ values obtained from a background 
subtraction are presented on Figs. 10-11 showing plots of projection angles Q x 
and Q z . 

The value of the parameter £ determined in the center of the volume (see 
Fig. 10), is in good agreement with estimations using Eq (5) giving £ esi m 23, 
t^sim ~ 22 (where £ est and £ S j m are the estimated and the simulated values 
respectively). The estimated values of £ for the tunnel are expected to lie 
within 2 — 5 because of the variation of parameter L along the tunnel. From 
the one dimensional plot (Fig 11a) one can see that the obtained results lie in 
the interval of the expected values. 

As one may expect, the summing of the signal bins along the tunnel will im- 
prove the signal to noise ratio (see Fig lib). So for thin heterogeneities with 
large lateral sizes it is not necessary to collect large statistics to become de- 
tectable (£ > 5 [27] suggested by formula (5)). So the approximate estimations 
for the statistics are good enough if the object's minimal lateral size and the 
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detector size are larger than (3 — 4) a v . 



7 Conclusions 

The muon intensity attenuation method to detect heterogeneities in matter 
volumes for the density length region of 5 — 50kg /cm 2 has been analyzed. 
Approximate analytical expressions for the signal to noise ratio and for the 
statistics collection time to achieve a given signal to noise ratio have been 
proposed. The energy integral intensity measurement is shown to be preferable 
than the differential energy intensity for the density length larger 5kg /cm 2 . 
The coordinate reconstruction and transfer uncertainties have been estimated 
using Monte Carlo simulations for different thicknesses of the absorber. They 
show a linear dependence on the distance from the detector and a smaller 
dependence on the muon low energy cuts. The minimum detectable cavity 
size has been shown to depend on the thickness of the matter volume along 
the muon track. 
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E (GeV) 



Fig. 1. Energy dependence of spectral index n{E): The solid line using parameteri- 
zation [15] and the dashed line using parameterization [14]. 
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Fig. 2. Energy dependence of ratios: F(n = constant) /F(n(E)) and 
&(n = constant) /&(n(E)) for three different zenith angles. The horizontal lines 
demonstrate 10% level of agreement. 
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Fig. 3. £ ratios (approximate/exact) dependent on L for different AL/L values. 




Fig. 4. The energy dependence of the power index m. The solid line is our fit based 
on CORSIKA simulated data at sea level, the dashed - the same at Mexico City- 
altitude and points, the experimental data for three different angles at sea level [21]. 
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Fig. 5. Available [2 1,22] large angle experimental zenith angular distributions for two 
representative energy values are shown as full circles and full squares, compared with 
the corresponding predictions of our fit (continuous and dashed, respectively) 
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Fig. 6. Experimental data[21] fit ratio for different zenith angles. The horizontal 
lines demonstrate 10% level of agreement. 
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7. Dependence of the <rc(full symbols) and ap (open symbols) on for different 
energy cuts of the detected muons. 
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Fig. 8. Simulation setup 
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Fig. 9. Distribution of detected muons as a function of projection angles considering 
that the elongated Smith tunnel and the hypothetical burial chamber are present. 
Layers of the pyramid and burial chamber (see text) are seen. 
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Fig. 10. Distribution of the sensitivity £ dependence on the projection angles. The 
tunnel is seen as a bright short strip on the top. The big chamber (see text) is seen 
very well. 
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Fig. 11. Distribution of the parameter £ in dependence on the projection angles a 
(a) and Q z (b). 
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